Discreteness-Induced Oscillatory Instabilities of Dark Solitons 



Magnus Johansson 1 and Yuri S. Kivshar 2 
1 Department of Physics and Measurement Technology, Linkoping University, S-581 83 

Linkoping, Sweden 

2 Optical Sciences Center, Research School of Physical Sciences and Engineering, The Australian 
National University, Canberra ACT 0200, Australia 



Abstract 

We reveal that even weak inherent discreteness of a nonlinear model can 

lead to instabilities of the localized modes it supports. We present the first 

example of an oscillatory instability of dark solitons, and analyse how it may 

occur for dark solitons of the discrete nonlinear Schrodinger and generalized 

Ablowitz-Ladik equations. 

PACS numbers: 03.40.Kf, 42.65.Tg, 63.20.Pw 
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Wave instabilities are probably the most remarkable nonlinear phenomena that may 
occur in nature 0. One of the first instabilities discovered for nonlinear models was the 
modulational instability, which is known to be an effective physical mechanism in fluids @ 
and optics for break-up of continuous modes into solitary waves. Also, the solitary waves 
themselves may become unstable, and the analysis of their instabilities is an important 
problem of nonlinear physics. Instabilities are known to occur for both bright M and dark 
H solitary waves of different nonintegrable nonlinear models. 

Recently, a new type of solitary-wave instability, oscillatory instability, has been found to 
occur for bright Bragg gap solitons in the generalized Thirring model M. Such an instability 
is characterized by complex eigenvalues, and its scenario is associated with a resonance be- 
tween the long-wavelength radiation and soliton internal modes which appear in the soliton 
spectrum when the model becomes nonintegrable |7j]. In spite of the fact that oscillatory 
instabilities appear often in dissipative models ||, their manifestation in continuous Hamil- 
tonian models is rare 0, and so far no example has been known for oscillatory instability of 
dark solitons. 

The aim of this Letter is twofold. First, we analyse what we believe to be the first 
examples of oscillatory instabilities of dark solitons, by considering the important cases 
of the discrete nonlinear Schrodinger (DNLS) and generalized Ablowitz-Ladik (AL-DNLS) 
models. We reveal two different scenarios for the dark-soliton oscillatory instability, which 
may occur due to either a resonance between radiation modes and the soliton internal mode, 
or a resonance between two soliton internal modes. Second, we demonstrate that even a weak 
inherent discreteness may drastically modify the dynamics of a nonlinear system leading to 
instabilities which have no analog in the continuum limit. 

First, we consider the well-known DNLS equation, 

ilp n + C(lp n+ i + 1p n -i) + \tp n \ 2 iJ n = 0, (1) 

where the dot stands for the derivative in time. Stationary localized solutions of Eq. ([I]) in 
the form ip n (t) = <p n e lM , where A = 2C cos k + (0(°)) 2 , may exist as dark-soliton modes with 

2 



the nonvanishing boundary conditions <p n — *> ±0(°)e lfcri (n — > ±00), provided the background 
wave is modulationally stable, i.e. for CcosA; < [TO. Without loss of generality, we can 
put the background intensity to unity, (p^ = 1. 

The structure of the dark-soliton modes of Eq. ([!]) has been discussed earlier flll,12 



Here, we consider the case C > 0, in which the background wave has k = it and is 'stag- 
gered', as shown in Figs. l(a,b). The transformation ip n — > (—l) n ip n immediately yields the 
corresponding 'unstaggered' modes (k = 0) for negative C. 

The dark-soliton modes presented in Figs. l(a,b) describe two types of stationary 'black' 
solitons m in a discrete lattice, the on-site mode (A-mode) centered with zero intensity at a 
lattice site, and the inter-site mode (B-mode) centered between two sites. These two modes 
can be uniquely followed from the continuous limit (C — > 00) to the 'anticontinuous' limit 
(C = 0). At C = 0, the A-mode takes the form n = (...- 1, +1, -1, 0, +1, -1, +1, ...), and 
it describes a single 'hole' in a background wave with constant amplitude and a 7r phase shift 
across the hole. Similarly, the B-mode takes the form <p n = (... — 1, +1, —1, —1, +1 — 1, ...), 
and it describes the lattice oscillation mode with a 7r phase shift between two neighboring 
sites and no hole. 

Linear stability of the A-mode for small enough C follows from Aubry's theorem (The- 
orem 9 in Ref. JL3|) relating the linear stability of multi-breathers to the extrema of 
the effective action which is a function of the relative phases a n of the single breathers. 
For the DNLS equation ([[]) with a positive nonlinearity, local minima of the effective 
action correspond to stable solutions. A perturbative expression of the effective action 
to order C 2 was obtained in Eq. (A10) in ||14|| . For small C, it is enough to con- 
sider the phase-interactions between nearest and next-nearest neighboring sites. Then, 
the lowest order contribution to the effective action from neighboring excited sites is 
X)n[2C|A| cos(a n+ i — a n ) — C 2 cos 2 (a n+ i — a n )}, while the contribution from the next-nearest- 
neighbor interaction is — C 2 J2 n ^n cos ( a n+i — «n-i) + 2C 2 cos(a no+ i — a no _i), where n is 
the site with the 'hole'. For the A-mode, ot n+ \ — a n = n, a n+ i — a n _i = 0, n 7^ no, 
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a„ 0+ i — a no -\ = n, so that the effective action has a local minimum and the mode is stable 
for small C. 

For the B-mode, there are two neighboring sites having the same phase at the center, 
which suggests that the corresponding extremum of the effective action is a saddle. Numer- 
ical investigation of the eigenvalue problem also shows that the B-mode is unstable for all 
C. 

Let us consider the linear stability of the dark-soliton modes for nonvanishing C. As- 
suming ip n {t) = [<f> n + e n (t)]e lAt , we obtain the linearized equations for e n 

ie n + C(e n+ i + e n _i) + 2|0 n | 2 e n + <p 2 n e* n - Ae n = 0. 

Writing e n = £ n + ii] n for real n yields 
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where for a system of iV sites M is a 2N x 2N matrix and H + and H~ are N x N matrices 
with time-independent coefficients, = [A — (2 =f l)(J)f}8ij — C(5 iy j + \ + (boundary 
conditions not explicitly taken into account). Linear stability is then equivalent to the matrix 
M having all its eigenvalues on the imaginary axis. 

When C — 0, all eigenvalues of M lie at zero except, for the A-mode, one complex 
conjugated pair at ±i corresponding to a mode localized on the 'hole'. When C is in- 
creased, the eigenvalues at zero spread on the imaginary axis creating a phonon band of 
extended states, which corresponds to a continuous spectrum for large N. Assuming e n = 
ae i(Kn-ut) + be -i( Kn -ut) yieldg the dispersion relation, u = ±^Jl6C 2 cos 4 (/t/2) + 8C cos 2 (/t/2), 
so that the eigenvalues corresponding to the continuous band lie between (at k — ir) and 
±W16C 2 + 8C (at k = 0). On the other hand, the pair of eigenvalues at ±i will move to- 
wards zero on the imaginary axis as C is increased. For small C, we can assume that the cor- 
responding eigenmode is almost completely localized at the hole, so that e no+ i ~ e no _i ~ 0. 
Since 0„ o = 0, we obtain that this mode oscillates with frequency A = 1 — 2C, so that the 
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corresponding eigenvalues are ±z(l — 2C), to the lowest order in C. Then, equating this 
result with the expression above for the edge of the continuous band gives an estimate to the 
value of C where these two eigenvalues coincide: C = l/y/3 — 1/2 « 0.07735. This agrees 
quite well with the exact, numerically obtained value C = C CT ~ 0.07647. At C — C cr , 
a Hopf-type bifurcation occurs, as two complex conjugated pairs of eigenvalues leave the 
imaginary axis and go out in the complex plane (see Fig. 2). Thus, an oscillatory instability 
occurs for the A- mode when C > C cr , with the instability growth rate A given by the real 
part of the unstable eigenvalue (see Fig. 2). For the B-mode, all eigenvalues lie at zero when 
C — 0. As soon as C is increased, one pair go out on the real axis and stay there for all 
C > 0. Thus, the B-mode is always unstable. 

For C > C cr , the imaginary part of the unstable eigenvalue moves towards zero (i.e, the 
oscillation frequency decreases) as C is increased. The real part of the eigenvalue vs. C 
is shown in Fig. 2 for two different system sizes, and compared with estimations of the 
maximum growth rate from direct numerical integration of the DNLS model for chains large 
enough to eliminate the influence of the boundaries. In all cases the instability is largest 
for C ~ 0.32 where A ~ 0.120. For larger C, there appear stable 'windows' for the finite 
systems, where the mode is stabilized by the boundaries. The location of these windows 
depend critically on the size of the system and the boundary conditions (in Fig. 2, periodic 
boundary conditions ipN+i — ipi are used). 

Similar 're-entrant instabilities' have been found for discrete breathers of the finite- 
width Klein-Gordon chains [15]. Qualitatively, the explanation is that for small systems, 



the eigenvalues corresponding to the 'continuous' spectrum are rather sparsely distributed, 
and the corresponding eigenvectors are localized over the size of the system. Thus, as 
the eigenvalue corresponding to the unstable mode approaches the imaginary axis when 
increasing C, it may find a 'hole' in the spectrum and join the imaginary axis for a while. 
Then, a further increase of C causes a collision with the next imaginary eigenvalue, and 
a new instability occurs. This procedure repeats itself until the eigenvalue has passed the 
phonon band eigenvalue which is closest to zero. After this the state is stable for all larger 



values of C (e.g., in Fig. 2 the state is stable for all C > 1.11 when N = 61, and for all 
C > 1.38 when N = 121). Increasing the system size implies that the eigenvalue distribution 
on the imaginary axis will be more dense, so that the stable windows will be smaller and 
completely disappear in the limit of the infinite system, where the unstable eigenvalue never 
joins the imaginary axis. 

For an infinite DNLS chain, the instability growth rate decreases in an exponential-like 
way to zero for large C, and thus indicates that the dark mode is unstable for all C > C cr . 
A quite good asymptotic fit is obtained with a stretched exponential, A ~ exp(— 6C 7 ), with 
7 ~ 0.7 (see inset in Fig. 2). We also investigated the case of varying (p^ to keep the 
complementary norm || constant when varying C . In that case, we found the asymptotic 
decay of the instability growth rate with C to be faster than purely exponential. 

Figure 3 shows the time evolution for two different values of C when the initial state is a 
slightly perturbed dark A-mode. A small perturbation has been chosen to be approximately 
in the direction of the unstable eigenvector. If instead a random perturbation is chosen, 
the qualitative behaviour is the same, except that there will be an initial small amount of 
radiation before the unstable internal eigenmode gets excited. The unstable eigenvector is 
always spatially symmetric around the central site, and therefore antisymmetric w.r.t the 
dark mode itself, since the latter is antisymmetric. Thus, the main effect of the instability 
is a transition from a black (zero intensity at the middle) to grey (nonzero intensity at the 
middle) soliton, as can be clearly identified at least for C > 0.3. This is similar to the 
instability scenario of dark solitons in some continuous models ||, except the oscillatory 
dynamics. The direction of the resulting grey soliton is determined by the sign of the 
perturbation projected on the unstable eigenmode. For large C, the grey soliton is almost 
black, moves slowly, and the radiation is small, while for smaller C the minimum value 
of the soliton intensity increases as well as its velocity and the amount of radiation. For 
small and intermediate values of the coupling parameter C, the resulting grey soliton decays 
continuously into radiation (see Fig. 3). 

It is important to study the oscillatory instability of dark solitons for other types of non- 
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linear lattices. Here, we consider the AL-DNLS equation iip n + C (ip n+ i +if) n -i) + fJ>\ip n \ 2,l Pn + 
\(n — l)|V , n| 2 ('0n+i + i>n-i) — 0, where < \i < 1. The case \i = corresponds to the 
integrable AL model, whereas the case \x = 1, to the DNLS model analyzed above. 

With the general form ip n (t) = (p n e tAt , and boundary conditions <p n — > ±0(°)e ifcn (n — > 
±oo), the frequency A is determined by A = 2C cos A; + (0 ( - o - ) ) 2 [/i + ( / u— 1) cos A;]. In particular, 
for the black mode (k = n) the relation becomes A = — 2C + ((p^) 2 , just as for the DNLS 
model. As above, we put (f)^ = 1 without loss of generality. 

In the AL-DNLS model, there is a lower limit of C for the existence of dark solitons, due 
to the instability of the background when the effective coupling changes sign (see, e.g., Ref. 
fl6||). This occurs when C + (// — l)(0^°^) 2 /2 = 0, so that, when (f>^ = 1, the dark modes 
exist only for C > (1 — fi)/2 (Fig. 4, dashed line). 

Considering ip n {t) = [<p n + e n (t)]e tAt as above, we obtain the linearized equations for 
the small perturbation e n , and the dispersion relation for the continuous spectrum (for the 
'staggered' mode with lim^i^oo 2 = 1): uj = ±2yV(/t)[l + u(k)], where z/(k) = (2C + ji — 
1) cos 2 (k/2). Therefore, the eigenvalues corresponding to the continuous band lie between 
(at k = 7r) and ±i2yJ(2C + fi- 1)(2C + fi) (at k = 0) for C > (1 - /i)/2. 

As above, an approximate expression for the onset of instability can be obtained by 
assuming that the internal mode is completely localized at the hole, so that e no+ i ~ e no _i ~ 
0. Then, the corresponding eigenvalues are ±i(l— 2C), to the lowest order in C. The collision 
with the phonon band edge is then expected to occur for C ~ (1 — 4/i)/6 + \/ fj, 2 + fi + 1/3. 
Comparing this result with the numerically obtained values of C for the bifurcation (see 
Fig. 4) shows that the instability for smaller [i occurs for smaller C than expected from 
this analytical estimate. The reason why the bifurcation occurs earlier is, that it in fact 
results not from a collision with the band edge phonon, but with a second localized mode 
that has bifurcated from the band edge slightly before (see the insets in Fig. 4). As a 
matter of fact, studying the bifurcation for the DNLS model very closely shows that also 
in this case the collision probably occurs with a second localized mode coming from the 
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phonon band. However, this localized mode only occurs extremely close to the bifurcation 
when fi is close to 1 (for /x = 1 we observe it for C > 0.076464, and the bifurcation occurs 
at C = 0.076468...). And since it will be very weakly localized, we are not able to say 
for sure whether it will be localized or not for the infinite system. But for smaller /i, this 
localized mode undoubtedly exists, and has the same spatial symmetry as the 'hole' mode. 
It describes an internal degree of freedom of the dark soliton, i.e. its internal mode. 

Otherwise, the qualitative scenario with stable and unstable windows for finite systems 
and instability all the way up to the continuum limit for the infinite system remains the 
same for all AL-DNLS lattices with < \i < 1. The exception is the limit of the AL model 
(/i = 0) which is known to be exactly integrable. At /i = 0, dark solitons are always stable 
for all C and their spectrum has a pair of eigenvalues exactly at zero corresponding to the 
exact translational invar iance. 

In conclusion, we have described, for the first time to our knowledge, the oscillatory in- 
stability of dark solitons. This new type of dark-soliton instability appears due to inherent 
discreteness of a nonlinear lattice model, and the universality of the instability scenario sug- 
gests that it should be also observed in other nonlinear models supporting dark solitons. We 
remark that the instability observed here may be regarded as an extension of the instabilities 
existing for small lattices [17] . 

M. J. acknowledges support from the Swedish Natural Science Research Council. Yu. 
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been reported at the CECAM Workshop on Nonlinear Localized Excitations in Condensed 
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FIGURES 



Fig. 1. (a),(b) Two types of 'staggered' dark-soliton modes in a lattice. Shown are the 
oscillation amplitude at each site (upper row) and the value \ipn\ 2 (lower row). 

Fig. 2. Real part of the unstable eigenvalues as a function of C in the model ([I]) for 
N = 61 (dotted) and iV = 121 (dashed) sites, and the growth rate (solid) for an infinite 
system, calculated for lattices of up to 15000 sites. Inset: asymptotic fit of the growth rate 
(points) with stretched exponential (line) (see text). Below: instability scenario shown as 
the evolution of the eigenvalues of M for iV = 121 sites. 

Fig. 3. Time-evolution of \ip n \ 2 with slightly perturbed dark A-modes as initial conditions 
for C = 0.30 (left) and C = 0.75 (right). Lower figures show the detailed dynamics for a 
few sites around the center. 

Fig. 4. Solid: Numerically obtained instability threshold C cr (fi) for the AL-DNLS model. 
Dashed: Minimum value of C for the existence of the dark localized modes. Insets: Examples 
of the eigenvalues at \i = 0.2, below (A) and above (B) the threshold curve C cr (/x). 
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